計算機で一様擬似乱数を生成する.
実数 \(a < b\) に対し,区間 \([a, b]\) 上の 一様分布 (uniform distribution) とは,確率密度関数 \[p(x) = \begin{cases} 0 & \text{if} \ \ x < a \\ \frac{1}{b - a} & \text{if} \ \ a \leq x \leq b \\ 0 & \text{if} \ \ x > b \end{cases}\tag{2.1}\] を持つ確率分布.\(\mathcal{U}[a, b]\) と書く.区間を指定せずに一様分布というと,\(\mathcal{U}[0, 1]\) を指す.
set.seed(1) # seedの設定
runif(3) # 長さ3の一様擬似乱数
## [1] 0.2655087 0.3721239 0.5728534
関数 runif の出力は一様乱数と認識しても問題ないか?
\(\rightsquigarrow\) スペクトル (Spectral) 検定,コルモゴロフ・スミルノフ (Kolmogorov-Smirnov) 検定 を適用.
# スペクトル検定
set.seed(1)
x <- runif(1000)
y <- runif(1000)
z <- runif(1000)
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)
立方体に均等に分布しているように見える.
# コルモゴロフ・スミルノフ検定
set.seed(1)
u <- runif(1000)
ks.test(u, punif)
##
## One-sample Kolmogorov-Smirnov test
##
## data: u
## D = 0.024366, p-value = 0.5928
## alternative hypothesis: two-sided
一様分布から乖離しているのであれば,\(p\)-値は小さくなるはずであるが,ここでは一様擬似乱数を観測として見たときの \(p\)-値は小さくない.
\(\rightsquigarrow\) 一様擬似乱数 (初期設定はMersenne-Twister法) は一様分布からの独立な乱数列と見ても顕著な差が見られない.
\(a, b, y_0, n\): 事前に決める整数.
\[ y_m = (a y_{m - 1} + b) \mod n \ (m = 1, 2, \ldots) \\ x_m = \frac{y_m}{n}\] として\(\{ x_m: m = 1, 2, \ldots \}\) を出力する.
u <- numeric(3*1e3)
y <- 1234
n <- 2^31 - 1
a <- 65539
b <- 0
for(i in 1:length(u)){
y <- (a*y + b) %% n
u[i] <- y/n
}
u_mat <- matrix(u, nrow=3)
x <- u_mat[1,]
y <- u_mat[2,]
z <- u_mat[3,]
plot_ly(x=x, y=y, z=z, type="scatter3d", mode="markers", size=0.1)
2次元平面に整列している様子が見て取れる.
以下,理論的な記述の際は乱数を扱い,実際の計算では擬似乱数を用いるが,擬似乱数も乱数と表記する.
ここで,\(\mathcal{U}[0, 1]\) からの一様乱数が生成できるとする.\(X \sim \mathcal{U}[0, 1]\) に対して,線形変換
\[ Y = b + (a - b) X \] を施すことで,\(\mathcal{U}[a, b]\) に従う乱数が生成できる.
\(\because\) \(\mathcal{U}[a, b]\) は累積分布関数 \[ F(x) = \frac{x - a}{b - a} \ (a \leq x \leq b) \] を持つが,線形変換してできた \(Y\) の従う分布も \[ \mathbb{P}(Y \leq x) = \mathbb{P}\left( X \leq \frac{x - b}{a - b} \right) = \frac{x - b}{a - b} = F(x) \ (a \leq x \leq b) \] となり,同じ累積分布関数を持つため.